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Using the variable phase method, we reformulate the Dirac equation governing the charge carriers 
in graphene into a nonlinear first-order diflerential equation from which we can treat both confined- 
state problems in electron waveguides and above-barrier scattering problems for arbitrary-shaped 
potential barriers and wells, decaying at large distances. We show that this method agrees with a 
^ ■ known analytic result for a hyperbolic secant potential and go on to investigate the nature of more 

experimentally realizable electron waveguides, showing that, when the Fermi energy is set at the 
. Dirac point, truly confined states are supported in pristine graphene. In contrast to exponentially- 

decaying potentials, we discover that the threshold potential strength at which the first confined 
state appears is vanishingly small for potentials decaying at large distances as a power law, but 
nonetheless further confined states are formed when the strength and spread of the potential reach 
^ ■ a certain threshold. 
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I. INTRODUCTION 



After being constructed in the 1920sA the variable phase method (VPM) was soon adapted to the theory of potential 
scattering by Morse and AUis- and has been expanded by numerous authors since as a method to solve scattering and 
' O ] confinement problems in quantum mechanics.'^"- The method has proven to be robust in the non-relativistic case and 
^ ■ has been applied to various physical problems, including an analysis of the statically screened Coulomb potentiali^ 
and the ionization of the electron-hole plasmaSiS in conventional semiconductors. The method has also been used to 
find scattering lengths of colliding atomsi^ and has been extended for use with non-local potentialsfiiii^ However, to 
the best of our knowledge, the method has never been applied to the two-dimensional modification of the Dirac- Weyl 
, equation describing the low-energy spectrum of charge carriers in graphene J^^— In this paper we study the presence 
^ • of bound modes within a smooth, confining electrostatic potential well in graphene in a waveguide geometry. 

° It is generally accepted that purely electrostatic confinement of charge carriers in graphene is not possible due 
■ to the effect of Klein tunnelingJ^"— At finite energy, carriers inside a potential well couple to states outside the 
] potential via same-energy states in the valence band. This is due to the conservation of chirality, which forbids 
. backscattering for normally incident particlesjii therefore the transmission probability is unity irrespective of the 
' barrier height. However, for particles possessing a finite longitudinal wavevector, which is the situation for particles 
T-H , that are not normally incident on the barrier, the effect is diminished. This suggests that a waveguide geometry, in 
^-H ■ which electrons with such a longitudinal wavevector propagate along the channel, may be appropriate for the creation 
^ " of conducting channels in graphene sheets. 

It has been previously demonstrated that the waveguide geometry for electrostatic potential wells does indeed 
lead to the lateral confinement of electrons in channels)^"— from which these confined modes provide an enhanced 
longitudinal conductivity. However, for energies away from the Dirac point the conductivity provided by the confined 
states is swamped by the sea of free carriers in the sample. This complicates the immediate application of graphene 
to digital electronics, which requires high on/off ratio current. It is therefore prudent to consider a graphene sheet 
close to the charge neutrality point, which can be easily achieved by adjusting the back-gate potential. 

Exact analytical solutions for a smoothly- varying electrostatic confinement potential are only known for a somewhat 
unrealistic hyperbolic secant potential at zero energy j^i^ and are also available at any energy for even less realistic 
square wells. The VPM allows one to probe the presence of confined states in any arbitrary potential channel as 
long as the potential decays faster than 1/x at large distances away from the central axis and is non-singular at finite 
distances. Therefore this method can act as a tool to determine the confined states in physically realizable electrostatic 
potentials. There is also continued interest in chiral tunneling and potential scattering problems in grapheneJ^i^ 
The VPM provides an efficient numerical tool for treating above-barrier reflection from arbitrary-shaped decaying 
one-dimensional potentials at non-normal incidence. 

The rest of this work is outlined as follows. After developing the formalism of the VPM in Sec. |TT] we show that 
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the method suggests a model top-gate structure should support zero-energy modes (Sec. IIII[) and captures a hitherto 
unknown major distinction between exponentially-decaying and more realistic power-law decaying potentials. Namely, 
there is a threshold potential strength for the appearance of a confined zero-energy state in an exponentially-decaying 
potential, whereas there is no threshold for a power-law decaying potential. This observation is confirmed by general 
analysis in Sec. lIVI We discuss the impact of our results in Sec.|Vl Appendices 1X1 and [Bl contain detailed derivations of 
the analytic results for two exponentially-decaying potentials, which were used to validate the proposed variable-phase 
method. 



II. FORMALISM 

In the low-energy approximation, charge carriers around a Dirac point in graphene are governed by the Dirac-Weyl 
Hamiltonian acting on a two-component wavefunctiou;^^ In the presence of an external electrostatic potential U{x) 
the Hamiltonian is given by^ 

H = VFiraxpx +<yyPy) + U{x) +TAaz, (1) 

where vf ~ c/300 is the Fermi velocity of the charge carriers, crx,y,z are the Pauli spin matrices and px,y — —iTidx.y are 
the linear momentum operator components. For completeness, and due to interest in the creation of confined states 
via band gaps in graphene, a mass term A is included in the final term in Eq. ([l}. Here r = ±1 denotes whether the 
low-energy approximation is around the K or K' points. Since the potential varies only in x, the momentum operator 
Py commutes with the Hamiltonian and so we may seek a solution in the form of the two-component wavefunction 
\['(x, y) = eKY){iqyy)['il)A{x),il>B{x)Y , where > is the longitudinal wave vector along the waveguide and 4'a{b) refer 
to the wave-function components associated with the inequivalent A{B) sublattices in graphene. These wave-function 
components satisfy the coupled first-order differential equations 

V{x) - T5)iipA, 

(2) 

V{x) + TS)iipB, 

where e = E/Tivf, E is the energy eigenvalue, V{x) = U{x)/hvF which we assume to be rapidly vanishing as a; — )■ ±oo 
and S = A/hvp. Notably, the sign of the potential is not important in considerations of zero-energy states in gapless 
graphene, as can be seen by reducing the system of Eqs. ([2]) to a second-order differential equation (see Appendix 1^ . 
Incidentally, there has been a recent rise in interest in zero-energy states in graphene, due to both a curiosity in edge 
states^^'^" and the possibility of observing Majorana zero-modes in this system. '^^■^^ 

Following the work of Babikov,— we write a wavefunction component as a superposition of transmitted and reflected 
waves 

^Pa{x) = C{x) (e"?.- + i?(x)e-*«^^) , (3) 

where C{x) and D(x) are known as the transmission functions and reflection functions respectively, and propagation 
is taken to be in the positive x direction. In addition to the introduction of the transmission and reflection functions, 
we will impose a condition on the first derivative of the upper wavefunction component, as we are free to do to 
completely define C{x) and D{x). 

Our choice of wavefunction in Eq. ([3]) is a modification of the free-particle solution which has the energy spectrum 
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e = ± (g^ + 9y + 1 which includes a band gap of 2J.~ The amplitude functions C{x) and D{x) have a natural 

interpretation if one considers the potential to possess cutoffs at a point xi to the left of the origin and X2 to the right 
of the origin. Then the barrier exists only in the region [xi, X2]', outside of this region the amplitude functions are not 
expected to be position dependent, thus C{x2) is simply the transmission amplitude and C{xi)D{xi) is the refiection 
amplitudeJ^ For a fast-decaying potential the amplitude functions tend to these position-independent values away 
from the center. Since the considered problem contains only a single electron and the electrostatic potential well, one 
may assert boundary conditions on the amplitude functions. The ingoing amplitude is unity, C{x —00) = I, and 
the refiection on the other side of the boundary is zero, D{x +00) — 0. The reflection coefficient^ can then be 
simply taken to be i? = \D{x — oo)|^. 

In implementing the VPM we make the following aforementioned ansatz on the first derivative of the upper wave- 
function component 

^i^Aix) = C{x) r Ae»9.- + Dix)^e-''^A . (4) 



dx 
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The ansatz indeed agrees with direct differentiation of the wavefunction in the hmit of large x, when the amphtudes 
become constants. However, this choice of ansatz is expedient because it always allows the transmission function 
to be eliminated later on (resulting in a differential equation purely in terms of the reflection function). Equating 
the right-hand side of Eq. Q with the full derivative of the wavefunction yields the following useful relation for the 
transmission function 

dx ^ ' dx e"J-^ + £'(a;)e-*9-^ " ^ ' 

Upon substituting Eqs. ^ and (HJ into the second of Eqs. ([2|) we find the lower component of the wavefunction to be 

^b{x) = ^^yl^^^^^g [irq. + iqy)e"^-'' + - Tg,)D(x)e-^«--] . (6) 

This is exactly the form for ij^B that one would expect from a flat potential profile, since it reflects the pseudospinor 
nature of the components. We may now employ the complete wavefunction and the first of the coupled Eqs. ([2]), and 
use Eq. (|4]) to eliminate the transmission function, resulting in a differential equation only containing the reflection 
function: 



dD{x) ^ fAjx) 
dx 2qx 



_dV{x) 
dx 



fB{x)+i[{e~V{x)f -e^) fA{x) 



(7) 



where the auxiliary functions fA.six) = iPa,b{x)/C{x). This is a nonlinear, first-order differential equation of the 
Riccati type and does not permit an analytical solution for an arbitrary potential. Numerical attempts to solve this 
equation will break down if e — V{x) + t6 = 0, which would occur for a potential in which the Fermi energy crosses 
both electron-like and hole-like regions. For real wavevectors, Eq. ([7]) describes scattering problems for Dirac fermions; 
reflection and transmission amplitudes can easily be found and the signature of resonant tunneling through the barrier 
is indicated by zeros in the reflection amplitude at certain values of energy and the parameters of the potential. For 
imaginary wavevectors, Eq. ([7]) allows one to tackle confinement problems for Dirac fermions, including analysis of 
bound and quasi-bound states, which are of special interest due to the phenomenon of Klein tunneling J^^— 

We are interested in probing the presence of bound modes in the potential at zero energy for the massless Dirac 
fermions of pristine graphene, thus we must set S — and = Hqy In the massless case the charge carriers exhibit 
valley symmetry, allowing us to confine ourselves to only one K point by setting t = 1 without loss of generality. In 
making the choice of complex wavevectors, the interpretation of the reflection probability from the potential barrier is 
that bound-state resonances are indicated by divergences (peaks) in R.^^ As a result, solving the differential equation 
numerically becomes almost intractable. These divergences are overcome by applying a convenient regularization of 
the form 

L>(.t) === tan[6'(a;)], (8) 

and rewriting the differential equation in terms of the phase function 9{x). Divergences in D{x) are mapped onto 
9{x) — ±7r/2, thus the numerical method becomes manageable. This is familiar from the Levinson theorem applied to 
the one-dimensional Dirac problem^^iiS^ in that the appearance of a new bound state corresponds to an accumulated 
scattering phase shift of tt. Taking = iqy, the first-order differential equation (O in terms of 9{x) for e = reduces 
to 



de 

dx 



-fix) 



1 dV{x) „„„,,,^,i_,„. , V{x 



l2 



cos[6'(a;)]e-'?«^ + -^-^f{x) 



V{x) dx 2qy 
fix) = cos[e{x)]e"^y'' + sm[e{x)]e''y' 



(9) 



This differential equation can be solved by employing the initial condition 9{x +oo) = 0, which arises due to 
zero refection amplitude in the region far behind the potential, and solving as an initial value problem. Eq. Q is a 
stiff differential equation so it requires implicit multi-step methods in order to accurately proceed to x — > — oo. The 
reflection coefficient R = \D{x — oo)P can then be foundi^ 



III. RESULTS 



We can now implement the VPM for massless Dirac fermions into physical problems. It is instructive to first 
demonstrate that this numerical method can reproduce analytical results from the literature. An exact result has 
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FIG. 2: (a) A schematic diagram of the experiment being considered. The top-gate is charged and produces an electrostatic 
potential in the graphene plane. The top-gate is treated as a wire of radius ro suspended above the graphene plane. The 
graphene layer is maintained at charge neutrality by charging the back-gate to Vtg. (b) Comparison of the realistic top-gate 
structure (dotted blue line) to the fitted hyperbolic secant model (solid red line), for h\/h2 — 0.9. 



been stated for an electron waveguide defined by the model potential V{x) — —Vo/cosh{x/d) in Ref. \2u\. However, 
since the derivation of this result"^^ has been left out of the widely available literature, we outline a procedure to 
obtain the exact solution in Appendix \X\ In this case, the condition for zero-energy confined states to be supported 
is iV'o'^l — Qyd + {n + i), where ti is a non-negative integer and the threshold for the appearance of the first mode 
is \Vod\ > i. For example, if we choose a potential with unit characteristic width {d = 1) and consider modes with 
propagating wavevector = ^ we expect zero-energy modes to appear at integer values of potential strength |Vo|. 
In Fig. m we plot the reflection coefficient R against potential strength for three values of qy , including qy — ^. The 
positions of the positive spikes in the reflection coefficient correspond to bound states, verifying the analytically-derived 
condition. 

The value of the proposed method is that it can be employed with any desired potential, including those not 
analytically solvable. Let us then consider a more physically realizable electrostatic potential, namely the potential 
created by a top-gate structure as in Fig. HJa). This structure can be formed by the deposition of a thin metallic 
strip onto an insulating layer of material on the graphene sheet. The insulating layer may be removed, for example. 
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by solvents, producing a so-called 'air bridge' ] ^^'^^ A model top-gate potential can be obtained, to first order, from 
the well-known method of image charges2S 



where hi and h2 are the distances between the metallic back gate and the graphene sheet and charged top gate, 
respectively, and T is matched to the top-gate voltage by 



where rp is the characteristic width of the top gate. For simplicity, we assume the dielectric constant of the dielectric 
layer between the graphene and metallic back gate to be equal to unity, which corresponds to the experimentally 
attainable case of suspended graphene.— The presence of a dielectric can be easily accounted for, again with the 
method of images;'**' however, it does not noticeably change the functional dependence of the potential, and especially 
its long-range asymptotics. To compare the top gate and hyperbolic secant waveguides we use a simple fitting 
procedure: fixing the potential maximum at the origin and the half width at half-maximumj^ This yields the following 
expressions for Vq, the absolute value of the potential at x = 0, and d in terms of top-gate parameters, for which the 
hyperbolic waveguide potential will best match that produced by the realistic top gate. 



For unit characteristic width {d = 1) of the hyperbolic waveguide we show an example of this fitting procedure in 
Fig. [21b), where we are using feasible experimental parameters;?^ hi « 2.7d, /12 « S.Od, and Vq « 2.94:eT /hvp. We can 
see that close to the origin the potentials fit well, since both functions are approximately parabolic. At large distances 
the behavior is very different, however; the secant hyperbolic waveguide is exponentially suppressed, whereas Eq. ()10|) 
falls as l/x^. This difference motivates this work because different behavior would be expected for states which 
possess quasiclassical turning points at significantly different distances from the origin. There is also a fundamental 
mathematical difference between the two confinement potentials under consideration; namely, the different behavior 
of the logarithmic derivative of the potential, which features in Ec^. In general the logarithmic derivative term 
will reach zero asymptotically for a power-law decay, whereas it will tend to a constant value for any exponential 
damping. Since an exponentially damped electrostatic potential is not realistic, we should be cautious when drawing 
general conclusions from such a model function. 

We compare the results of the analytic solution'^" and the more realistic potential in Fig.lSja). As expected, in both 
cases new propagating states appear fairly regularly with increasing potential depth. The appearance of bound states 
for the hyperbolic secant waveguide follows a linear relation between Vq and qy and new states appear equally spaced 
according to the analytical relationship. The top-gate also produces new states following a roughly linear relationship 
for large propagating wavevectors. Such a similarity is to be expected when one notes that an electron possessing a 
large propagating wavevector will have quasiclassical turning points close to the origin, where the potentials possess 
similar functional forms. 

The appearance of confined modes for the realistic top-gate potential does not follow a linear relationship for lower 
Qy] instead as Qy decreases, the values of Vq, at which the first states appear, fall off dramatically. In fact, contrary to 
the hyperbolic secant waveguide (and indeed the exponentially decaying cusp waveguide which also yields an exact 
solution, see Appendix |B]) , which possesses a finite threshold value, for the realistic top-gate waveguide case we find 
a complete absence of a threshold, as shown in Fig. ^h). This could be attributed to the logarithmic derivative 
term, which acts drastically different in the two cases (as we shall see in Sec. I IV I) such that only potentials with 
exponential tails have a threshold potential strength. This should lead to an enhanced conductivity in a graphene 
sheet for waveguide-type devices even for vanishingly weak potential channels. Interestingly, the zero threshold for a 
power-law decay of the confining potential does not apply universally to charge carriers in graphene. It was shown 
recently in the radial geometry that electrostatic potentials without exponentially decaying tails may indeed possess 
a threshold for the appearance of the first zero-energy state4ii^ 




(10) 
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(11) 




(12) 



IV. GENERAL ANALYSIS OF POWER-LAW DECAYING POTENTIALS 



It has been stated in Ref. [1^ (see also Appendix |X| and shown numerically in Sec. IIIII that for the hyperbolic 
secant waveguide there exists a threshold potential strength, below which there are no bound states. Here we show 
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FIG. 3: (Color online) The appearances of bound states for the top-gate potential and the fitted exactly-solvable hyperbolic 
secant waveguide, with h\/h2 = 0.9 and d — 1. Boxes and crosses indicate propagating modes for the model and realistic 
top-gate potential correspondingly. The left panel shows the large-scale plot for qy € [0,0.5]. The right panel contains a 
logarithmic plot for small qy for the realistic potential only, indicating that the threshold for the power-decaying top-gate 
potential is vanishingly small. 



that this property does not extend to potentials that decay as a power law, such as the considered model top-gate 
potential Eq. ([T0|. due to the logarithmic derivative of the potential tending towards zero asymptotically. 

First, let us consider the case of a general exponentially-decaying potential V{x) = Vq/ cosh^{x/d) with parameter 
t > 0. As the logarithmic derivative will always tend to a constant asymptotically, in the limit x — >■ ±oo, it follows 
from Eqs. ^ that 

and the equation for ipAix) is obtained upon replacing Qy with —qy in Eq. (I13p . Hence the normalizable solutions 
describing the first emergent state of both wavefunction components in the region a; > are 

i^Bix) - e-9«% tPAix) - e-(««+*/'')^ (14) 

and accordingly in the region a; < 

tPeix) - e^iy-*/'^>, ^^(a;) - e'^''^ (15) 

Thus, we arrive at the condition qyd > t to ensure a non-trivial normalizable solution in the limit of large negative 
X and so it is apparent in this scenario we do indeed have a threshold, as we expect from our numerical solutions in 
Sec, mil A particular case of the hyperbolic secant potential is considered in Appendix \X\ The presence of a threshold 
in the product of potential strength and spread for the appearance of the first confined mode can also be shown 
explicitly for another exponentially-decaying potential supporting analytic zero-energy solutions which is considered 
in Appendix [Bl 

Now let us turn our attention to the equation governing the behavior of the lower wavefunction component ips (x) 

in the potential V{x) = Vb/(1 + fr)^^^, where the decay is characterized by p > 1. In the limit x ±oo Eqs. ([2]) 
yield 

^i^Bix) + ^^^b{x) + - qA ^Bix) = 0, (16) 

dx^ X ax \ X "J 

and a similar equation for '\1)a{x) can be written by changing qy to —qy in Eq. (jl6p . Both of these equations can 
easily be reduced to a confluent hypergeometric differential equation, known as Kummer's differential equation4^ The 
square-integrable solutions describing the first emergent state for both wavefunction components when a; > are^ 

^B{x)r.e"iy\ x1^a{x) ^ e~'iy^ U{p,p,2qyx), (17) 

where lA denotes the confluent hypergeometric function of the second kind, or Tricomi function. The corresponding 
square-normalizable solutions for x < are'*'* 



i)B{x) ~ e'""^ U{p,p, -2qyx), 2Pa{x) - e^" 



(18) 
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The terms including the Tricomi functions decay quickly at large x (in the region considered) for all decay strengths 
p and propagating wave vectors qy, as can be seen via a series expansion at infinity.*- The functional forms of the 
wave- function components at large distances, Eqs. (|17 | I18 |) . indicate that no conditions are required to be imposed 
upon the longitudinal wave vector qy and hence there is no threshold of potential strength at which bound states first 
appear. 

V. SUMMARY 

The variable phase method is elegant because it allows one to solve the Schrodinger equation directly for physical 
quantities, such as the reflection coefficient or the scattering phase, rather than needing to extract these properties 
from the wave function.— We have shown that this method can be extended to the Dirac-Weyl equation. We expect 
this method to become broadly used for modeling top-gate devices in the considered waveguide geometry, as well as 
for calculating barrier transparency at non-normal incidence, by computing reflection and transmission coefficients. 
The VPM also offers significant practical advantages: whilst it is always possible to reduce two coupled first-order 
differential equations into a decoupled second-order equation, the VPM allows one to decouple into a single first-order 
equation, from which a useful physical property can be immediately obtained. Of course, the nonlinear nature of the 
equation is of little importance when integrating numerically. 

In conclusion, we have modified the variable-phase method so that it can be used to study quasi-one-dimensional 
problems, such as electron channels, for the quasi-relativistic charge carriers of graphene. The method has been 
validated by reproducing the exact result for a hyperbolic secant waveguide. We have gone on to find, for a model top- 
gate electrostatic potential with experimentally obtainable parameters, that there is no potential strength threshold 
for bound state emergence, so the first confined mode exists for an arbitrarily weak power-decaying potential. With 
increasing potential strength further confined states appear in steps of applied top-gate voltage of the order of tens 
of millivolts. 

In fact, we have found that while there is always a threshold potential strength for the appearance of the first 
confined zero-energy state in channels defined by a potential with exponential tails, there is no such threshold in 
potentials decaying as a power law. This suggests that both the square waveguide model and indeed models with 
exponential tails are somewhat unsatisfactory, because they miss an important physical property necessarily present 
in more realistic top-gate defined potentials. The presence of zero-energy modes for all realistic power-law decaying 
waveguide potentials can be related to the well-known problem of non-vanishing conductivity in graphene for the case 
of smooth disorder. Our results provide an additional argument in favor of the so-called resistor network model for 
minimum conductivity.— 

Electrostatically-defined waveguides are free from apparent shortcomings of graphene nanoribbons such as strong 
back-scattering produced by edges and technological difficulties in controlling nanoribbon parameters. Moreover a 
top-gate controlled potential prevents edge-related scattering and can serve as an analog of an electrostatically-defined 
carbon nanotube with an enhanced mean free path. We expect truly confined states to be observable in charge neutral 
pristine graphene without the need of magnetic fields, strain engineering or the introduction of band-gaps, a potential 
way forward in the pursuit of graphene devices aimed at digital manipulation of current. 
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Appendix A: Exact solution of the hyperbolic secant waveguide 

This appendix illustrates a method of solving the massless Dirac-Weyl equation for the case of a hyperbolic secant 
electron waveguide at zero-energy, following previous work in Ref. [36| . The wavefunctions and threshold conditions for 
fully-confined states in this potential were stated in Ref. 20] without derivation. Since the aforementioned threshold 
conditions were used as the major benchmark for checking validity of the formulated VPM, we found it appropriate 
to provide here the analytic solution of the hyperbolic secant problem. 
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Upon setting E = and decoupling Eqs. ([2]) into a second-order differential equation for a single wavefunction 
component ipBix) only, we obtain 



(f , , 1 dV{x) d , , , ( , ,-y Qv dV(x) , , 



(Al) 



Considering the one-dimensional potential V{x) = — Vq/ cosh(a;/c?) and parameters gj, > 0, d > and making the 
change of variable ^ — tan'h{x/d) yields 



(A2) 



where w — \VQd\ and S — qyd. This is a known differential equation which should be reduced using the following form 
of the solution^ 



1 



(A3) 



where the function is to be found and p and q are subject to the conditions 

iqiq -l) + 2q + a + b + c = 0, {p - q)[2{p + q) - 1] = c, 

a = —oj'^, b = S, c = ui'^ — S'^ . 

Satisfying the conditions for p and q and making a further change of variable k = (^ + l)/2 leads to a differential 
equation in the well-known Gauss hypergeometric form 



d^ 

k{k - 1)^?7(k) -f 



(2p + 2q-l)K- [ 2p 



1 



dn 



r]{K)+ [ip + qf+a] v{k) = 0, 



from which one can write down the unnormalized wavefunction solution as 

= (1 + 0^(1 -0*2^^1 (p + q + u;,p + q-uj;2p+l^ 



(A4) 



(A5) 



where p = ^^^^ ~ j-, Q = ^^3^ + \ and n is a non-negative integer. To avoid a singularity at ^ = ±1 one obtains the 
condition that uj — n > ^. This puts an upper limit on n, such that for a channel of given parameters there may exist 
only a finite number of distinct propagating states. Termination of the hypergeometric series via \Vod\ — qyd — n + ^ 
gives the condition for which confined states are supported. The upper wavefunction component, tJja, can be obtained 
from the coupled Eqs. ^ and has a similar form to ipB- 



Appendix B: Exact solution of the exponentially decaying cusp waveguide 



This appendix details the solution of the massless Dirac equation at zero-energy for an exponentially decaying cusp 
electrostatic potential^^ defined by V{x) = Vb exp(— ja^j/ci), with parameters qy > 0, d > 0. The aim of this appendix 
is to confirm our general statement on the presence of threshold in the product of potential strength and its spatial 
extent for the appearance of the first confined state in exponentially-decaying potential. Besides, to our knowledge, a 
solution of the two-dimensional Dirac- Weyl equation for this potential has not been reported in the literature before. 

The second-order differential equation to be solved for wavefunction component tpsix) follows from Eq. (|A1[) 



dx^ 



yjBix) + 8gn{x)~Mx) + (l^o'e-'l^'l/'^ + sgn(x) | - q^) 4>b{x) = 0, 



(Bl) 



where sgn(a::) is the signum function of the coordinate x. This equation is of a standard form and yields the following 
solution in terms of Bessel functions of the first kind"*^ 



(B2) 
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where the first term in the solution is present only for a; > and the second term only for a; < 0. Again, the upper 
wavefunction component tpj^ can be found via the coupled Eqs. Continuity of both wavefunction components at 
X = requires Ci — C2 and yields the constraint 

J,,d^i/2 {\Vod\) = ±J,„<j+i/2 {\Vod\) . (B3) 

Equation (|B3[) is a transcendental equation, and graphical or numerical solutions display a threshold value of \Vod\ = 
0.785... for which the potential can support its first bound state, and indeed the values of potential strength at which 
further bound modes appear. 
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